Combining quantum and classical density functional theory for ion-electron mixtures. 
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We combine techniques from quantum and from classical density functional theory (DFT) to 
describe electron-ion mixtures. For homogeneous systems, we show how to calculate ion-ion and ion- 
electron correlation functions within Chihara's quantum hypernetted chain approximation, which 
we derive within a DFT formulation. We also sketch out how to apply the DFT formulation to 
inhomogeneous electron-ion mixtures, and use this to study the electron distribution at the liquid- 
solid interface of Al. PACS numbers:71. 22. -f 1,61.10.-1,61.20. Gy,61. 12. Bt 



I. INTRODUCTION 

Density functional theory (DFT) has proven itself a remarkably successful tool in condensed matter physics [Q. 
The foundations were laid in 1964, when Hohenberg and Kohn (HK) Q proved that the ground state energy of any 
quantum mechanical system could be described as a functional of the one-body density only. Subsequently, Kohn 
and Sham [|j developed an orbital based method which could be applied to electronic systems, and Mermin 0| 
extended the HK proof to finite temperatures, opening up the possibility of using DFT to calculate the free-energy of 
a statistical mechanical system. Since then many different practical methods have been developed to apply DFT to 
electronic problems, with countless applications in condensed matter physics, chemistry, and biology g. Examples 
of electronic structure techniques relevant to this paper are the ah initio molecular dynamics (AIMD) method of Car 
and Parrinello |^ and the orbital free ah initio molecular dynamics (OF- AIMD) scheme of Madden and co-workers 
0. In a parallel development, the finite temperature version of DFT has been widely used to study classical systems, 
with myriad applications to phase-transitions and the theory of fluids (in its broadest sense), see e.g. |^,^ for reviews. 

However, the classical and quantum versions of DFT have rarely been combined - the two fields have largely 
developed independently. In this paper we will attempt to bring together these two strands of DFT into one unified 
formalism, with the aim of applying them to liquid metals. These are ideal systems for which to test a mixed 
quantum-classical DFT approach since they can be viewed as binary mixtures of classical ions and quantum electrons. 

To begin, we first define a two-component free-energy functional 

F[pi,p,] = Ff[pi] -f Fl%,] -f F^^[pi,p,], (1) 

split into the usual way between the ideal free energies of the ions, Fj'^, and the electrons, Fl'^ , and the excess free 
energy F^^ , which is a functional of both one-body densities pi{v) and /Oe(r) . This unique intrinsic free-energy 
functional functional obeys the variational principle; 



5F[pi,pe 
6 Pair) 



+ '/'a(r) = 0, (2) 



where the (/>Q(r) are the external potentials which uniquely induce the one-body densities Pa{^)- This principle is 
subject to the constraint that the integral of Pa{v) = Na, the total number of particles of species a — {e, /}. We will 
apply this free-energy functional both to homogeneous and to inhomogeneous liquid metals. 

In section || we use the DFT formulation to derive the quantum Ornstein Zernike(QOZ) equations. From these 
we can derive a set of integral equations, called the quantum hypernetted chain approximation (QHNC) by Chihara 
, which self-consistently solves for the ion-ion and ion-electron pair correlation functions in a homogeneous liquid 
metal. The great simplification of the QHNC approach is that the original many-centre electronic problem is reduced 
to an effective one-centre problem, which is much easier to solve. As a specific application, we study the ion-electron 
radial distribution function for a set of simple metals. 
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In section III we extend this DFT formulation to inhomogeneous fluids. By using a simple version of this functional, 
we compare the electron density profile at the liquid-solid interface of Al to the recent simulations using OF-AIMD 
101. We end with some conclusions. 



II. APPLYING TWO-COMPONENT DFT TO HOMOGENEOUS LIQUID METALS 

A. Quantum Ornstein Zernike equations 

DFT provides a nicely unified route to the correlation functions of homogeneous fluids. To begin, we flrst show 
how to derive the QOZ relations by defining the direct correlation function as the second functional derivative of the 
excess free-energy defined in Eq. (|]): 

where (xa^a)"^ is the inverse susceptibility matrix of the full system, and (x^^)^^ is the inverse susceptibility matrix 
of the ideal system. These relationships between direct correlation functions and inverse susceptibilities are called 
the QOZ equations, and they hold for an arbitrary number of components. If any components are classical, then the 
susceptibilities can be connected to the pair-correlations through the fluctuation-dissipation theorem : 

limxa/3(fc,0) = -/3(p°p°)i/25„0(A:), (4) 

where the Sapik) are the structure factors, defined in the usual way [|l3|;0, and we have taken the homogeneous 
limit and the Palpha^ are the homogeneous limits of the densities PQ(r). Specializing to the homogeneous limit 
for two components, with one component being classical (the case at hand), the QOZ equations reduce to a set of 
relationships between the direct correlation functions and the the ion-ion and ion-electron pair correlation functions as 
well as the electron-electron susceptibilities. The latter are still quantum-mechanical, and cannot be simply connected 
to electron-electron pair correlations. 



B. Quantum Hypernetted Chain Approximation 

To make further progress, one needs a way of solving the QOZ relations for a liquid metal. We sketch the derivation 
here, for more details we refer to reference jl^] and to the original work of Chihara [ p^Jl7| . As a first step, we use 
a quantum version of the Percus trick [^8| to relate the homogeneous two-body pair-correlation functions to the one- 
body inhomogeneous density around a single classical particle fixed at the origin. For the ion-electron pair-correlation 
function we fix an ion at the origin to find: 

g^r)^^, (5) 

Pe 

where pe(r\I) is the (interacting) valence electron density A similar relationship can be found for the ion- ion 

pair-correlation function, but this trick cannot be used for the electron-electron pair-correlation function, since one 
cannot "fix" an electron at one position. The next step is to solve for the interacting one-body electron-density. To do 
this, one can use the Kohn-Sham trick ||], namely that there exists a single-particle external potential v°^{r) which 
will induce in a non-interacting system the same one-particle density p{r) as is found in the full interacting system. 
This external effective potential, felt by the non-interacting electrons or non-interacting ions, follows from the Euler 
equations. If one further makes a functional Taylor expansion around the equilibrium homogeneous densities, then 
the following effective ion-ion and ion-electron potentials result: 

vfi{r)=vMr)^^Y.Py J Co.^i\r-r'\)h^i{r)dr' + ^B^jir), (6) 

where the Ca-yir) are the homogeneous limits of the direct correlation functions defined in Eq. (|^), and the Percus 
trick was used to rewrite (p^(r|w^/) — p'^) in terms of the correlation functions h^jir) = g-yi(r) — 1. The remaining 
higher order terms are lumped into the so-called bridge functions Ba0{r), well known in the theory of liquids po| , |l4| . 
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Again, we note that no such potential can be derived for the electron -electron correlations. With that caveat in mind, 
our formulation, derived from a DFT approach, is still in principle exact. 

To make progress, however, some approximations must be made, in particular to treat the electron-electron corre- 
lations. We follow Chihara, who made a series of such approximations to derive the QHNC equations |0. The most 
important one is: 

1. The valence electron correlations are treated in the jellium approximation. I.e. the effect of the ion-ion and ion- 
electron correlations on the electron-electron correlations are ignored, and the electron-electron direct-correlation 
function is written as: 

Cee(fc) = -l3Vee{k)[l " G^f (7) 



where is the local field factor for the electron gas This has the great advantage that the QOZ 

relations (y) now reduce to only two coupled equations for the ion-ion and the ion-electron correlations. While 
this is the key approximation that makes the QOZ relations tractable, it is also the most important and 
uncontrolled approximation in the QHNC scheme. The fact that the effective pair potentials derived from a 
linear response scheme are very sensitive to the exact form of the local field factor also testifies to the 

importance of the approximations in Eq. to the physics of a liquid metal. 

QHNC also uses an implicit separation of the valence electrons from the core electrons. This approximation 
is closely related to the jellium approximation discussed above, and is only good when the core and valence 
electron energy levels are well separated so that the ionic correlations don't have a significant effect on the core 
states; it is likely to break down near resonances. 

The other approximations used in QHNC are much better understood, and are often encountered in either their 
electronic structure or liquid state theory context . In descending order of importance they are (For a 
more detailed discussion see ref. [p^.): 

2. The local density approximation (LDA) is used for the one-centre ion-electron problem. This approximation 
is well known and understood in electronic-structure theory ||2l|] . Although in principle this approximation is 
similar to approximation 1 for Gee(Q), since the exchange and correlation energy of an inhomogeneous electron 
gas are approximated by the local energy density of jellium, there are a number of reasons to believe the LDA 
is rather accurate ||^. Various improvements to the LDA exist, but these are not thought to be very important 
for the simple metals we study. 

3. The ion-electron bridge function _B/e(r) is set to 0, which is related to the hypernetted-chain approximation, 
well known in liquid state theory Jlif . 

4. The ion-ion bridge function Bn (r) is approximated by the bridge-function of a hard-sphere reference state. This 
is called the RHNC or MHNC approximation in liquid state theory |ljj2^, and is generally very accurate. 

5. The bare ion-ion potential is taken to be purely Coulombic. Again, this is an approximation whose limitations 
are well understood . 

In summary then, the QHNC approximation achieves the following very important simplification for the electronic 
problem: the original many-centre electronic problem has been reduced to an effective one-centre problem by replacing 
the direct effect of the ions with an effective external potential, given by Eq. (|^), that depends self-consistently on 
the ion-ion correlations. This is depicted schematically in Fig. 0. The main advantages are that the ion-ion and 
ion-electron correlations emerge naturally and on the same footing and that the calculations are much more rapid 
than full ah initio simulations. By using DFT to derive the QHNC equations, their origin and meaning becomes more 
transparent. For details on the (non-trivial) numerical implementation of the QHNC we refer to refs. |lq-0|. 
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FIG. 1. Schematic picture of the QHNC approximation; The original many-centre problem is 
reduced to an effective one-centre problem. 



C. Applications of the QHNC to ion-electron correlations 



The ion-electron radial distribution function, defined as the conditional probability of finding a valence electron a 
distance r away, given that there is an ion at the origin, can be written as: 



where n{r) is the so-called "pseudo-atom" density, which, when superimposed according the ion-ion radial-distribution 
function gii{r), gives the correct value of the total electron-density. The ion-electron radial distribution function and 
the related pseudo-atom density, calculated with the QHNC approach, are shown in Fig. ^ for a series of simple 
metals. At short distances the probability of finding an electron a distance r away is equal to the pseudo-atom 
density, but further away the effects of the other surrounding pseudo-atoms kick in. Note that the QHNC is an 
all-electron calculation; the oscillations near the core are correctly reproduced in contrast to the traditional 
AIMD approaches, which rely on pseudopotentials. The ion-electron correlation functions calculated with the QHNC 
compare very well for the cases where AIMD simulations are available. However, it turns out that for most of the 
metals in our set a much simple linear response formalism also performs remarkably well |23[ | . The surprising accuracy 
of the linear response arises from a quantum interference effect between two length-scales: the Fermi wave- vector and 
the core-size of the ions, which makes the non-linear response terms much less important than one might naively 
expect p^]. The accuracy of the QHNC also benefits from this effect. 




(8) 
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FIG. 2. The ion-electron radial distribution functions as obtained from the QHNC approximation 
(sohd lines), OF-AIMD |34[||l (open circles) and Car-Parrinello AIMD [|§ (crosses). The dashed 
lines represent the pseudo-atom density n{r)/ p^. 



The only metal in our set where the QHNC has difhculty is Ga. We found earlier |15| that an ad- hoc change in the 
local field factor Gee{q) seems to improve matters. We also found that the d-electrons are close to a resonance. These 
two points indicate that approximation 1, discussed in the previous section, begins to break down for Ga. 



III. APPLYING TWO-COMPONENT DFT TO INHOMOGENEOUS FLUIDS 

In this section we explore the possibilities of studying inhomogeneous metals within a two-component DFT formu- 
lation. The advantage of treating metals directly as as electron-ion mixtures in this way is that we bypass the need 
for effective (densitydependent) ion-ion potentials. Modern classical DFT methods use input from the homogeneous 
liquid state the QHNC is ideally suited to provide such input for a two-component DFT. 

One of the simplest approximations for the DFT of an inhomogeneous system, first proposed by Ramakrishnan and 
Yussouf (RY) [^,0 , corresponds to truncating at second order a functional Taylor expansion of the excess free-energy 
around the homogeneous phase: 

AF-[p,,pe] = -EE /^^ / dr'Ap^{r)CMplPp;\r-r'\)App{r'), (9) 

a p •' 

where AF'^^ = F'^'^[{pa{r)}] — Fq^{{p^}), and the direct correlation functions Ca^iPa: P°p\ |r — r'|) are those of the 
homogeneous liquid phase. For metals, the QHNC provides Cii{r) and Cie{r), while Cee{r) is fixed by the jellium 
approximation. Together, these correlation functions completely determine the excess free-energy of Eq.(^. 

To obtain the total free-energy, one then adds the ideal contributions of the electrons and the ions. For the ions, 
this ideal free energy functional has the form: 
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PF^Pi] = / dvpiiv) {In [p/(r)A3] - 1} , (10) 

but for the electrons the exact form is not known, and various approximations must be made. To preserve the 
advantages of the current DFT approach, one needs a kinetic energy functional. A number of forms have been 
proposed; some are remarkably accurate, see e.g. and references therein. These kinetic energy functionals also vary 
in the ease with which they can be implemented in a full two-component DFT calculation Pql . 

Classical DFT has frequently been applied to the freezing transition and the fluid-solid interface of simple fluids 
1^ . The two-component DFT described above provides a basis upon which to build similar applications for metallic 
systems. 

As a first application, we attempted to use the two-component RY formalism described above to calculate the 
freezing transition of Al. A very similar approach was already implemented by one of us [ p9| to study the freezing 
transition of H. However there a simpler semi-empirical prescription was used for Ciejr). Since H does not have a 
core radius, the cancellation of higher order response terms found for simple metals [E3[ does not hold. This implies 
that the correlation functions and also the freezing transition of H should be more difficult to treat accurately than 
would be the case for simple metals. Other properties of liquid 77, such as the metal-insulator transition are also 
much harder to treat. Partially for these reasons, we would initially expect that the simple metals are actually better 
systems on which to test rudimentary DFT's. 

The one-body densities of the solid phase are parameterized by sums over the reciprocal lattice vectors (r.l.v.) {G}: 

Pa(r) = p°^^^exp[zGT]. (11) 

{G} 

For the ions a further, well established, simplification is to take the real-space density pi{v) as a sum of Gaussians, 
centred on the lattice sites, which implies S^f = exp [—0^/4^] . With this density parameterization and approximation, 
the full free-energy functional for the difference between the solid and the liquid state free energy, taken from Eqs. ^ g 
and |l^, reduces to: 



IvT 



{G} 



I E 'ierP%e{\G\) ^ E '(?ep"eCle{\G\). (12) 
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{G} {G} 

Here the Cq,/3(|G|) are the Fourier transforms of the direct correlation functions for the liquid state, and the primes 
mean that the term G = should be omitted in the sum. The first term between brackets is the ideal free-energy of 
the ions (^0|), written out in terms of the Gaussian parameter 

This free-energy can then be minimized through the variational principle (^) to solve for the free-energy difference 
between a given solid crystal structure and the fluid phase. However, when we attempted this for Al, with the direct 
correlation functions given by the QHNC, we found that the free energy of the solid did not have a stable minimum. 
This is a common problem encountered when applying a DFT to the freezing transition (see e.g. ref ||3l[]): If the 
full direct correlation functions are used, then it is very hard to stabilize solid phase. If instead one uses the DFT 
for a hard-sphere reference system, adding the attractive interactions as a perturbation, very good agreement with 
experiment and simulation is recovered. Examples include phase diagram |Q and the fluid-solid interface of the 
LJ fluid and a one-component treatment of the freezing of Al p^ l, based on an effective pair potential. This 
suggests that a successful 2-component DFT treatment of the freezing transition also needs to pass via this route. 
Unfortunately this quickly reduces to a theory similar to the one-component DFT, which means introducing the 
effective pair potentials we are trying to avoid in the first place. 

Another application of the two-component DFT formalism is to examine the fluid-solid interface of Al, for which 
Besson and Madden |l^ have recently completed extensive OF-AIMD simulations. They found non-trivial behaviour, 
with the ion-ion density profiles decaying from an ordered structure to a smooth liquid structure across several atomic 
layers. Similarly, they found that the valence electron density also showed oscillations which decayed into the bulk 
liquid phase. In principle, a two-component DFT would be ideally suited to study this problem, since the ionic and 
electronic density profiles come out on equal footing. Unfortunately, the difficulties encountered when trying to treat 
the freezing transition also make studying this problem very difficult. Instead, we tackle a simpler problem: Given 
the ionic density profile, can we use our DFT to calculate the electron density profile? 

If the ionic density profiles are given, which fixes the ^f, then the electronic which satisfy the variational 
principle (^) follow from Eq. (p^): 
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G 



Cf5/e(G) 



(13) 



(/3/xL(G))+Cee(G))- 

where we have made an additional approximation, also made in [29j, namely that the ideal electron contribution 
AFg'^[pe] is expanded to second order in (linear response). This is consistent with the spirit of the RY DFT 
approach, and greatly simplifies the minimization of the electronic degrees of freedom, since they are now linearly 
coupled to the ionic ones. Note that taking the lowest order approximation to the ion-electron direct correlation 
function, Cie(r) w —(3vie{r) in Eq. ( p^ would be equivalent to a standard linear response calculation. Since we 
use the full C/e(r), non- linear response terms are included, although one would need to use the full kinetic energy 
functional to be completely consistent with the RY DFT approach of Eq. (n2|). The results of using Eq. (O) are 
shown in Fig. || as a function of the distance z along the interface. (See ref. for details of parameters)"^ The 
QHNC results are in fact less accurate than a simpler linear-response approximation which replaces c/e(G), with 
an empty-core Ashcroft pseudopotential At first this may seem surprising, but a careful look at the inset 

of Fig. y shows that the difference results from the fact that c/e(G) is considerably smaller at G ~ 1.8, where the 
strongest ionic scattering occurs. This difference stems in part from the fact that the true pseudopotential in Al is 
known to be highly non-local ||2l[] because of the p character of the bonds, suggesting that a local formulation, such 
as the simple linear response DFT we employed, will not work well together with the QHNC calculation which makes 
no local pseudopotential assumptions. Non-linear response effects are also expected to be more important at the 
highly inhomogeneous fluid-fluid interface. In contrast to the QHNC, the pseudo-potential contains an empirically 
adjustable parameter which effectively includes non-local and non-linear effects. To achieve accurate results with a 
DFT based on the QHNC, one needs to treat the two approaches at comparable levels of approximation. As is often 
encountered in condensed-matter physics, mixing different levels of approximation does not work very well. Two ways 
of improving this would be (1) to use QHNC with the improved local pseudopotentials used in the OFMD or 
(2) to use a more sophisticated DFT approach. Although far from complete or completely satisfactory, this first, and 
rather primitive, example of a two-component DFT used for a liquid-solid interface suggests that our DFT approach 
could in principle be fruitfully used to study fluid-solid interfaces in simple metals. 



0.005 




-0.0025 

40 

z (au) 

FIG. 3. Electron density profiles from the OF-AIMD simulations (solid lines), compared to 
the QHNC-DFT results (dashed lines) given by Eq. (p^), and linear response theory (dotted lines) 
with an empty-core pseudopotential with Rc — l.lSSao- Inset, comparison of cie{k) (solid line) and 
the pseudopotential vi^{k) (dashed line), both in Hartrees. 
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IV. CONCLUSIONS 



In conclusion, we have shown how DFT provides a unified formahsm from which to derive the QOZ equations, and 
the related QHNC approach, first pioneered by Chihara. The QHNC reduces the many-centre electronic problem to 
an eflFective one-centre one. This dramatically reduces the computational time needed to calculate ion-electron and 
ion-ion correlation functions. It also has the advantage that no explicit use of effective pair potentials or pseudo- 
potentials is needed. The most important approximation in the QHNC approach is to treat the electron-electron 
correlations as those of jellium. It is not yet clear under which conditions this approximation begins to break down. 
With this caveat in mind, the other approximations entering the QHNC are reasonably well understood. We used 
the QHNC to calculate the ion-electron radial distribution functions for a number of simple metals, finding very good 
agreement with ab-initio molecular dynamics calculations where these are available. 

We have also discussed an exploratory application of the simple RY DFT approach to liquid metals as two-component 
electron-ion mixtures. We found that describing the freezing transition suffers from a similar problem to that found 
for classical simple fluids: the solid phase does not develop a stable minimum if the full direct correlation functions 
are used as input. This is rather disappointing. Slightly more (partial) success was found when we attempted to 
use the DFT to describe the inhomogeneous electron density at a liquid-solid Al interface. Although this result was 
only a partial application of the DFT, since the ionic density profile came from the simulations, it does suggest that 
developing a full fledged DFT could be very fruitful. Finally, it is clear that our rudimentary DFT approach is not yet 
accurate enough, and could be improved in a number of ways. We are attempting to generalize more sophisticated 
DFT schemes like the MDWA psf or GELA § to the two-component ion-electron problem. We are also studying 
the effect of using different local held factors and different kinetic energy functionals. 
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